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ROBUST  IMPLEMENTATION  OF  LEMKE'S  METHOD  TOR  THE 


LINEAR  COMPLEMENTARITY  PROBLEM 


J.A.  Tomlin 


Abstrac 


This  note  discusses  techniques  for  implementing  Lemke’s 
algorithm  for  the  linear  complementarity  problem  in  a numerically 
robust  way  as  well  as  a method  for  recovering  from  loss  of  feasi- 
bility or  singularity  of  the  basis.  This  recovery  method  is  valid 
for  both  positive  semi -definite  M matrices  and  those  with  positive 
principal  minors.  It  also  allows  a user  to  start  from  an  advanced 


basis  for  such  problems 


/ 


l 

’ 
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Introduction 

While  Lemke's  algorithm  [L]  for  the  linear  complementarity 
problem  (LCP)  is  well  established  in  the  theoretical  sense  and  to 
some  extent  as  a computational  procedure,  it  has  not  come  into 
universal  use  for  solving  LCP's  in  the  same  way  that  the  simplex 
method  has  for  linear  programs.  This  is  due  to  several  reasons. 

To  begin  with,  some  LCP's  with  special  structure  can  be  solved  more 
efficiently  by  other  special  purpose  algorithms  (see  [1]).  It  is 
also  apparent  that  some  problems  (such  as  quadratic  programs  with  a 
relatively  small  proportion  of  quadratic  variables),  which  can  be 
set  up  as  LCP's  and  solved  by  Lemke's  method,  might  be  more  efficiently 
solved  by  direct  optimization  methods.  The  question  of  applicability 
of  the  algorithm  will  not  be  considered  here  however.  The  problem 
we  shall  be  concerned  with  is  the  stability  of  Lemke’s  algorithm  in 
the  face  of  numerical  error  and  its  ability  to  recover  in  the  face 
of  these  errors. 

The  revised  simplex  method  (see  e.g.,  [5]),  when  it  breaks 
down,  usually  does  so  in  one  of  two  ways  — a loss  of  feasibility 
in  the  solution,  or  the  current  "basis"  is  detected  to  be  singular 
(or  nearly  so)  when  an  attempt  is  made  to  invert  it.  Both  of  these 
situations  are  quite  easily  remedied.  In  the  first  case  the  algorithm 
can  revert  to  Phase  I of  the  simplex  method  to  restore  feasibility. 

In  the  second  case  dependent  columns  can  be  removed  from  the  basis 

and  replaced  by  logical  (unit)  vectors  to  give  a new  non-singular 
basis.  The  resulting  basic  solution  will  now  almost  certainly  be 
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"Revised"  Form  of  Lemke's  Algorithm 

Familiarity  with  both  the  Lemke  algorithm  and  the  revised 
simplex  method  is  assumed,  but  we  make  the  following  observations 
for  the  sake  of  clarity. 

The  linear  complementarity  problem  is  customarily  expressed 
in  the  form  (see  e.g.,  Cottle  and  Dantzig  [2]): 

Find  w and  z such  that 


w = q + Mz 


(1) 


T 

w,  z > 0 , wz=0 

A new  non-negative  vector  e which  is  positive  in  each 
component  corresponding  to  a negative  q^,  and  a new  variable  Zq, 
are  added  to  the  problem,  which  can  be  stored  internally  in  the 
form: 


Iw  - Mz  - ezQ  = q . (2) 

Since  column  e is  immediately  introduced  into  the  basis  to  provide 
a non-negative  basic  solution,  and  remains  there  until  termination, 
it  is  important  that  e only  contain  the  minimum  number  of  elements 
required,  otherwise  the  basis  factorization  or  inverse  representation 
may  become  extremely  dense. 
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infeasible  and  Phase  I must  be  carried  out  to  restore  feasibility. 

It  is  clear  that  the  ability  of  the  simplex  method  to  recover  in 
this  way  is  a result  of  its  two  phase  nature  and  the  admissibility 
of  any  basis  as  a starting  point  for  recovery.  These  features  also 
enable  the  user  to  start  the  simplex  method  from  any  purported  basis. 

The  Lemke  algorithm  is  also  vulnerable  to  loss  of  feasibility 
and  singularity  of  the  basis  (we  assume  a revised  simplex-like 
implementation).  However  since  there  is  no  "Phase  I"  and  the  basis 
is  required  to  be  of  a special  nature  (an  almost  complementary  set 
of  columns)  it  is  not  immediately  clear  how  these  misfortunes  can 
be  overcome.  Similarly  we  cannot,  in  general,  start  with  a given 
basis  unless  it  happens  to  be  feasible  and  almost  complementary. 

These  are  severe  restrictions  for  practical  computation  for  all  but 
the  smallest  and  simplest  problems.  In  view  of  these  difficulties, 
three  new  algorithms  for  TCP’s  were  developed  by  Rarick  [6]  which  did 
allow  for  a two-phase  procedure  and  which  successfully  solved  problems 
with  hundreds  of  variables.  At  present,  only  two  of  these  algorithms 
have  been  proved  to  converge  and  then  only  for  semi -definite  matrices. 
Since  Lemke's  algorithm  is  known  to  solve  (theoretically)  a wider 
class  of  problems,  it  seems  worthwhile  to  try  to  find  robust  procedures 
for  implementing  it  and  to  devise  methods  for  dealing  with  infeasibility 
and  basis  singularity. 


I^t  B be  the  current  almost  complementary  feasible  basis 
for  (2)  at  any  iteration  and  suppose  (wp,  z ) is  the  complementary 
pair  of  variables  of  which  one  has  just  become  non-basic.  Ihe  steps 
of  an  iteration  are  then  extremely  simple: 


! 


! 


!•  If  Zq  has  become  non-basic  or  < e,  for  some  small  €, 
terminate. 


2.  Determine  which  of  w^  or  z is  to  enter  the  basis  (the 

P p 

complement  of  the  variable  which  has  just  left  the  basis)  and  compute. 


Qt  = B~\i  or  a = -B-1m 
P P 


th 


where  m is  the  p column  of  M and  u the  p^1  unit  vector. 


3.  Compute 


Sc  . ft 

= mi“ 1 
v £^>01 


where  0 = B~  q 


i ! 

H 

- • 

i ' 


Record  the  new  non-basic  pair  and  update  0 and  the  representation 


of  B . (Either  the  standard  product  form  or  seme  decomposition 
of  B) . Go  to  1. 


The  basis  is  of  course  periodically  reinverted  or  refactorized. 
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Stable  Implementation 

Unlike  the  simplex  method,  the  variable  to  enter  the  basis  is 
uniquely  determined  by  the  previous  iteration.  There  is,  in  general, 
no  possibility  of  rejecting  the  column  because  of  unsatisfactory 
pivot  size  and  using  some  other  column.  Because  of  this  we  must  be 
as  careful  as  possible  in  choosing  the  pivot  in  step  3 of  each  itera- 
tion. The  first,  and  most  obvious  precaution  is  to  use  Harris's 
pivot  row  selection  technique  [3]  to  choose  the  largest  available 
pivot  in  the  event  that  there  is  a choice  of  pivots  which  maintain 
feasibility  within  a given  tolerance  e . Briefly  restated,  this 
involves  a two  stage  process: 


3a) 


cp  = 


a. 

1 


mm 

>*P 


Pi 


+ c 


0 


ct. 

1 


3b)  Choose  Pr/ar  such  that 


a 

r 


maxfa.l  a.  > e . 
1 1 - P* 


< cp) 

a.  - 

i 


where  is  an  absolute  pivot  tolerance. 

As  a further  check  we  may  impose  a relative  pivot  tolerance 

test  on  a . Let  a = maxi  a.|  and  €_  be  a tolerance.  Then 
r max  1 i1  R 

if  ar  > the  pivot  is  accepted.  If  not,  the  choice  is 

suspended  and  the  basis  is  reinverted  so  that  a will  be  recomputed 
as  accurately  as  possible  (or  practical) . If  the  recomputed 
again  fails  the  relative  pivot  test  it  must  ..uw  be  accepted,  but 
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with  greater  confidence.  Clearly  if  is  too  small  the  test  is 

ineffective  and  if  it  is  too  large  reinversion  will  he  carried  out 

-8 

more  frequently  than  necessary.  We  have  found  a value  of  lCf  to 

-15 

be  satisfactory,  when  the  relative  machine  precision  is  about  10 

Much  of  the  concern  about  accuracy  can  be  removed  by  using  one 
of  the  stable  techniques  for  updating  triangular  factors  of  the 
basis  now  available,  for  example  Saunders'  efficient  version  [7]  of 
the  Bartels-Golub  method.  We  should  also  point  out  that  use  of  a 
triangular  updating  technique  is  particularly  beneficial  on  sparsity 
grounds,  since  the  e column  in  (2)  persists  in  the  basis  and  may 
be  dense.  In  this  case  we  expect  it  to  be  assigned  as  one  of  the 
right-most  columns  of  U in  the  LU  decomposition  of  B produced 
by  almost  any  sparsity-preserving  technique.  In  this  position  it 
will  have  minimal  effect  on  the  build  up  of  non-zeros  during  the 
updating  process. 

As  a final  measure  we  may  periodically  check  the  relative 
error  in  3 by  performing  one  iteration  of  iterative  refinement 
with  the  current  representation  of  B \ If  this  relative  error  is 
unsatisfactory  reinversion  may  be  called  for. 

Matrix  Scaling 

If  the  matrix  M is  badly  scaled  or  ill-conditioned  the 
possibility  of  numerical  difficulties  will  be  magnified.  It  there- 
fore seems  worthwhile  to  have  scaling  facilities  available.  Again 
the  situation  here  is  not  as  straight-forward  as  it  is  for  the  simplex 
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method,  as  it  is  not  clear  that  all  matrices  M and  right  hand 
sides  q for  which  the  LC'P  is  solvable  by  Lemke's  algorithm  preserve 
this  property  when  the  problem  is  scaled.  Fortunately  this  solv- 
ability is  preserved  for  two  large  classes  of  matrices  — those  which 
are  positive  semi-definite  and  matrices  with  positive  principal 
minors.  If  we  let  R = diagCr.^  . ..,  r^)  and  C = diagfc^,  ...,  c^) 
be  the  diagonal  matrices  of  positive  row  and  column  scales,  then  it 
is  clear  that  if 


M'  = RMC 


then  M’  retains  semi -definiteness  and/or  positive  principal  minors 
if  M has  these  properties.  Hence,  the  scaled  problem 

w'  = q1  + M'z ’ 

w',  z'  >0  (w,)^zr  = 0 

where  w ' = Rw,  q'  = Rq  and  z'  = C~^z  remains  solvable. 

There  remains  some  ambiguity  however,  about  the  applicability 
of  scaling  to  the  more  esoteric  classes  of  matrices  for  which  Lemke's 
algorithm  is  known  to  be  effective.  In  particular,  it  is  not  clear 
whether  the  additional  column  e must  be  row  scaled  in  the  same  way 
as  the  matrix,  or  constructed  after  scaling.  Further  research  on 
this  question  is  called  for. 


The  choice  of  scaling  methods  for  M is  essentially  the  same 
as  for  linear  programming.  An  earlier  study  of  LP  scaling  [8] 
indicates  that  the  method  due  to  Curtis  and  Ried  is  very  effective 
and  we  recommend  it  here. 

Recovery  Procedure 

Despite  all  precautions  there  are  bound  to  be  occasions  when 
the  algorithm  breaks  down,  leading  to  an  infeasible  solution  or  a 
singular  "basis."  Rather  more  frequently,  a user  may  wish  to  supply 
a basis  which  turns  out  to  be  infeasible  or  even  singular  but  is 
close  in  some  sense  to  the  solution.  Both  of  these  difficulties  can 
be  overcome,  at  least  for  the  two  large  classes  of  M with  which  we 
are  principally  concerned  — positive  semi -definite  and  positive 
principal  minors. 

As  soon  as  either  error  condition  is  detected,  the  e column 
is  dropped  from  the  basis  and  replaced  by  the  unit  column  u^  corres- 
ponding to  the  Wp  variable  in  the  non-basic  complementary  pair. 

If  the  resulting  basis  remains  singular,  dependent  columns  m.  are 

J 

dropped  from  the  basis  and  replaced  by  their  complementary  unit 

columns  u..  This  process  must  lead  to  a non-singular  complementaiy 
J 

basis  — that  is  one  containing  a (possibly  empty)  set  of  columns 
-Mj  and  a complementary  set  of  columns  of  the  unit  matrix  Ij. 

Let  B be  the  resulting  basis. 


Now  consider  the  transformed  problem 


w = q + Mz 


w,  z ' >0  w £ = 0 


where  w = B ‘'w,  q = B ^q  M=B^M  and  w,  £ are  a complementary 
permutation  of  w,  z. 

By  the  definition  of  B,  M is  a principal  transform  of  M 
and  if  M is  positive  semi -definite  and/or  has  positive  principal 
minors  so  does  M (see  Cottle  and  Dantzig  [2]  Theorems  8 and  9). 

Thus,  if  the  original  problem  (1)  is  solvable  by  Lemke's  algorithm, 
so  is  the  transformed  problem  (5) . All  that  remains  is  to  construct 
a new  column  e with  positive  elements  corresponding  to  the  negative 
elements  of  q.  This  is  then  assigned  a new  variable  £^  and  added 
to  the  problem,  the  pivot  row  (and  new  nonbasic  pair)  being  chosen 
to  give  a new  almost  complementary  feasible  solution. 

It  is  important  to  note  that  we  continue  in  practice  to  deal 
with  the  original  M and  q and  B not  M and  q,  hence  we  require 
e also  in  terms  of  the  unit  basis.  This  is  of  course  available  as 
e'  = Bi,  and  e is  overwritten  with  e'  in  the  interval  representa- 
tion of  the  problem  (2) . 
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To  illustrate  the  fallibility  of  this  recovery  procedure 
consider  the  small  example: 


Constructing  e = (1,  1)  and  applying  Lemke's  algorithm  leads  in 
two  iterations  to  the  correct  solution  w^  = 9/11,  wQ  = 0,  z^  = 0, 
z2  = 5/11.  If  however,  we  pivot  on  to  produce  the  principal 
transform: 


the  resulting  M matrix  is  the  negative  of  a Minkowski  matrix  and 
the  algorithm  terminates  after  the  first  iteration  in  the  almost 
complementary  ray  zQ  = 2 + e,  = 1 + 56,  z2  = 6,  w = w?  = 0. 

'//hat  recovery  procedures  (if  any)  are  applicable  to  problems 


which  lie  outside  the  classes  we  have  considered  therefore  remains 
a topic  for  further  research. 


: _ 


Conclusion 

With  the  numerical  techniques  and  the  recovery  procedure 
discussed  in  this  note  there  is  no  reason  that  Lemke's  method  should 
not  be  used  for  a wide  class  of  large  scale  linear  complementarity 
problems.  With  the  current  exception  of  LU  updating  we  have  imple- 
mented these  techniques  in  our  code  LCPL  [9]  and  found  it  highly 
satisfactory.  So  far  our  experience  has  only  been  with  problems 
of  up  to  a little  more  than  500  variables;  however,  some  of  these 
had  led  to  difficulties  of  the  type  described  with  other  implementa- 
tions. 
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